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The purpose of this paper is to present a multilayer primitive equations model for ocean 
dynamics in which the velocity and buoyancy fields within each layer are not only allowed 
to vary arbitrarily with horizontal position and time, but also with depth — linearly at 
most. The model is a generalization of Ripa's inhomogeneous one-layer model to an 
arbitrary number of layers. Unlike models with homogeneous layers, the present model is 
able to represent thermodynamics processes. Unlike models with slab layers, i.e. those in 
which the layer velocity and buoyancy fields are depth-independent, the present model 
can represent explicitly the thermal- wind balance within each layer which dominates at 
low frequency. In the absence of external forcing and dissipation, energy, volume, mass, 
and buoyancy variance constrain the dynamics; conservation of total zonal momentum 
requires in addition the usual zonal symmetry of the topography and horizontal domain. 
The model further possesses a singular Hamiltonian structure. Unlike the single-layer 
counterpart, however, no steady solution has been possible to prove formally (or Arnold) 
stable using the above invariants. It is shown here that a model with only two layers 
provides an excellent representation of the exact gravest baroclinic mode phase speed. 
This suggests that configurations with only a small number of layers will be needed to 
tackle a large variety of problems with enough realism. 



1. Introduction 

Approximate inhomogeneous layer primitive equations models in which the velocity 
and buoyancy field s are allowed to vary only in the horizontal position and time ( Dronkcrg 
ll969HLavoill972l) have been very extensively exploited in ocean modeling. For instance, 
a model with only one "slab" active layer (i.e. floating on top of a denser homogeneous 
quiescent infinitely deep layer) has been shown useful to study processes near the ocean 
surface whi ch, in addition to wind drag, is subiect to nonuniform heat and freshwater 
fluxe s (e.g. lY oimg"l994^ 'Young Ch^ll99,4 lFuka,ma,chi ^7^1199,4 iRinalll 994 iR.piedl 
ll997tlScott fc Wi Umott 2002). The main motivation for the use of these type of models 
is the observation that in the upper part of the ocean a vertically well-mixed layer sepa- 
rated from the underlying heavier fluid by a well-defined pycnocline often develops. Also, 
the slab models are numerically very efficient, particularly in the study of meso- and 
submesoscale ocean dynamics as they can resolve, for a given comp utational cost, finer 
horizontal scales than fully three-dimensional models llEldevikll200'l . However, applica- 
tions are not restricted to the study upper-ocean processes only. In these ap plicat i ons, a 
stack of homogeneous and inhomogeneous active layers is usually considered: iRipal l)l993() 
has derived a model with an arbitrary number of active slab layers. Examples of these ap- 
plications include one-dimensional (in the horizontal) modeling efforts, typically oriented 
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to provide insight into the basi c physics underlying the heat and salt balances in semien- 
closed seas (e.g. Rim. ,199 1: Beron- Vera fc Ripa,.,2002: Ripa,2001.): and two-dimensional 
regiorial ocean modeling efforts fe.g. ISzo eke fc RichmanI 119841 And^rson&^dcCreanJ 
_19856-.'McCrearv fc Kundu"l988': McCreary et aL"l989'. 'l99l'. 'l993"l 996QHBeiejll997r 
Beier fc Ripa .1999: R0cd fc Shi 1999; Palacios-Hernandez et al. 200^V Other applica- 
tions include the study of biophysical int eractions (e.g. iMcCreary et oil Il9966l I2OOII: 
Hood et aL'2003'). equatorial dynamics fe.g. 'Anderson fc McCreary'l985a':'McCrear v fc Yul 
1992|) . and ENSO predictability fe.g. .Schopf fc Cane 1983: Balmascda e t al, 1994). Fur- 
thermore, in the widely used Miami Isopyc nic-Coordinate Oc ean Model (MICOM) the 
upper layer is chosen to be a slab-like layer iBl eck et QiJll992l) . 

Despite their widespread use, the slab models are known to haye several limitations 
and deficiencies iRioal l| 19991) : (a) they cannot represent explicitly the thermal- wind bal- 
ance which dominates at low frequencies; and (6) they have a zero-frequency mode in 
which the layer thickness and buoyancy changes without changing the flow, and whose 
unlimited grow th cannot be p revented by the invariants that constrain the (inviscid un- 
forced) system. lEldevikI l|2002|) has recently pointed out also that the slab models cannot 
satisfactorily represent frontogenesis, a process for which the thermal- wind balance is 
fundamental. 

To cure the slab model limitations and deficiencies, Ripa (1995, hereinafter referred to 
as R95) proposed an improved closure to partially incorporate thermodynamic processes 
in a one-layer model. In addition to allowing arbitrary velocity and buoyancy variations 
in horizontal position and time, Ripa's model allows the velocity and buoyancy fields 
to vary linearly with depth. Ripa's model enjoys a number of properties which make it 
very promising. For instance, (a) it represents explicitly the thermal- wind balance at low 
frequencies; and (6) with only one layer the free waves supported by the model (Poincare, 
Rossby, midlatitude coastal Kelvin, equatorial, etc.) are a very good approximation to 
the first and second vertical modes in the fully three-dimensional model. 

In this work I present a generalization of Ripa's model to an arbitrary number of 
layers, including two possible (mathematically equivalent) vertical configurations (§13). 
The goal is to generalize Ripa's model, allowing enough flexibility to treat more com- 
plicated problems than those that can be tackled with only one layer. Several aspects 
of the generalized model are discussed in f|31 These include: its conservation laws; the 
Kelvin circulation theorems derived from the model; its Hamiltonian structure; the lack 
of formal stability theorems; results on vertical normal modes and remarks on baroclinic 
instability; and the incorporation of forcing in the model equations. Section ^closes the 
paper with some concluding remarks. 



2. The generalized model 

Consider a stack of n active fluid layers with thickness hi{'x.,t)^ i = 1, • • • ,n, where x 
is the horizontal position and t the time (figure^. The geometry can be either planar 
or spherical; in the former case the vertical coordinate, z, is perpendicular to the plane, 
whereas in the latter it is radial. The total thickness is h{x.,t) = J^j hj{x,t). The stack 
of inhomogcneous layers can be either limited from below by a rigid bottom, z — /io(x), 
or from above by a rigid lid, z — —ho{x). The usual choice in the rigid lid case is 
ho = 0; however, laboratory experiments are often design to have a nonhorizontal top 
lid. The remaining boundary in the rigid-bottom (resp., rigid- lid) configuration is a soft 
interface with a passive infinitely thick layer of lighter (resp., denser) homogeneous fluid 
of density p„+i. Although vacuum {p^+i = 0) is the typical setting in the rigid-bottom 
configuration, the choice p^+i 7^ is useful in the study of deep flows over topography. 



Pn+1 = const. 

U„ + l = 0, pn+l = 
hn + 1 OO 




Pn+1 = const. 

U„+l = 0, Pn+1 = 

h„+i OO 
(a) (b) 



Figure 1. Sketch of the layered modeL The two possible vertical configurations of the model 
are rigid bottom (a) and rigid lid (b). Within each layer the velocity and buoyancy fields not 
only vary arbitrarily with the horizontal position and time, but also linearly with depth. 



Following R95 closely, I consider the ansatz 

Ui(x, a, t) = Ui(x, t) + cru^(x, i), (2.1a) 
i?»(x,a,i) =i9»(x,t)+aC(x,i), (2.1fe) 

for the ith-layer horizontal velocity and buoyancy fields. Here, 

n + l 



da a, (2.2) 

1 



which is the vertical average of any variable a(x, tr, t) within the ith layer. The ith-layer 
buoyancy is defined as 

i3i(x, cr, t) := ±g [p,(x, a, t) - p„+i] /p,, (2.3) 

where the upper (resp., lower) sign corresponds to the rigid-bottom (resp., rigid- lid). Here, 
g is gravity, Pi{x, a, t) = Pi{x, t) + cpf (x, t) is the (variable) density in the «th layer, and 
Pr denotes the (constant) reference density used in the Boussinesq approximation. A key 
element for the derivation of the n-layer version of Ripa's model is the definition of the 
rescaled vertical coordinate a so as to vary linearly from ±1, at the base, to =Fl, at the 
top, of the ith layer (figure l^j): 

±z=: /i,_i(x,t)-hi^/i,(x,t), (2.4) 

where 

i 

/i,(x,i) :=/io(x)+^/ij(x,i). (2.5) 
[Henceforth an upper (resp., lower) sign will correspond to the rigid-bottom (resp., rigid- 



Figure 2. Vertical coordinate choice. Within each layer the rescaled vertical coordinate a varies 
linearly from ±1, at the base, to =pl, at the top. The upper (resp., lower) sign corresponds to 
the rigid-bottom (resp., rigid-lid) configuration of figure 1. 



lid) configuration.] Physically admissible buoyancy values, i.e. everywhere positive and 
monotonically increasing (resp., decreasing) with depth in the rigid-bottom (resp., rig-lid) 
case, are such that 

>C>0, 1?, ^C+C+i- (2.6) 
Notice that -dl — ^nfhi, where nf{x,t) > is the square of the instantaneous Brunt- 
Vaisala frequency within the ith layer. 

In order to obtain the equations for the n-layer version of Ripa's model one must 
proceed as follows. 

(a) Substitute the ansatz p. 1(1 in the (inviscid unforced) fully three-dimensional, i.e. 
exact, primitive equations defined in ho < ±z < ho + h with vertical velocity w = 
±(9t -I- u • V)(/io + h) a.t z ^ ±(/io + h) and w — ±u • V/iq at z = ±ho, and kinematic 
pressure p = at z — ±(/io + h); here V denotes the horizontal gradient. 

(b) Replace all occurrences of by its vertical average (i.e. ct^ i-^ |) to preserve the 
linear vertical structure within each layer. 

(c) Collect terms in powers of a and equate them to zero afterwards. 
The resulting equations are: 



= 0, (2.7a) 
{Di^r = 0, (2.7b) 
dth, + V • h,u, = 0, (2.7c) 



D,u, + /z X u, + Vp, = 0, (2.7d) 
(D,u,)" + /z X < + (Vp,;)" - 0. (2.7e) 
Here, / is the Coriolis parameter; z is the vertical unit vector; 

B;^={dt + • V)a + Ihr^V ■ /i,<a", (2.7f) 
(D.a)'^ = {dt + u, ■ V)a" + < • Va, (2.7g) 

are the mean and a components of the material derivative of any vector a(x, cr,t) = 
a(x, t) + cra'^(x, t) in the zth layer; and 

n 

J =1+1 

{Vp,Y = + \hVd, + <V/i,_i, (2.7i) 

which are the mean and a components of the ith-layer pressure gradient force. 

System ()2.7|l consists of In evolution equations in the In independent fields {hi^-di, 
f?^, Ui, uf , • • • , /i„, i?„, 1?^, u„, u^). The couphng among different layer quantities is pro- 
vided by the last terms on the right hand side of the pressure forces (|2.7h .i). It is impor- 
tant to notice that the dynamics in both the rigid-bottom and rigid-lid configurations is 
described by system ((2.7|l : no double signs are necessary. The latter must be taken into 
account, however, in the computation of the total pressure in the ith layer, which, up to 
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the addition of an irrelevant constant, is given by Pj.pi ± p„_(.]^gz, where 

n 

K = i (1 + a) h,d, - i(l - a^)h^d^ + hjdj. (2.7j) 

The multilayer slab model follows from (|2.7|) upon neglecting uf and -d^ . Ignoring uf 
in (|2.7|l results in a model with uf = but ^ which provides a generalization of 
Schopf & Cane's (1983) intermediate layer model. Alternatively, omission oidl in system 
H2.7|l gives a model with 'd^=0 but uf ^ that can be thought of, assuming a rigid- 
lid configuration, as the generalized multilayer nonsubinertial mixed-layer approximation 
of lYoung l|l994) . A initial state with uniform buoyancy inside each layer {-di = const, 
and = 0) and vanishing velocity vertical shear (uf = 0) is preserved by H2.7|l : as 
a consequence, unlike the above reduced models, the homogeneous layers model follows 
from H2.7|l as a particular case, just as it does it from the fully three-dimensional model. It 
is important to notice that the homogeneous layers model is exact for a stepwise density 
stratification; however, it is not able to accommodate thermodynamic processes which 
are particularly important in the upper part of ocean. 



3. Discussion of several aspects of the model 

3.1. Conservation laws 

In a closed domain, D, conservation of the jth-layer volume, mass, and buoyancy variance 
is enforced, respectively, because of (12.7b '). 

dt{hA) + V • ih,m, + iz9>n = 0, (3.1) 

and 

dt (/i^^) + V-h, (^u, + p.i?>r) = 0. (3.2) 
The total energy (sum of the energies in each layer) is also preserved in a closed domain 

as 

3 3 

where 

E, := i/i,u2 -I- \h,{ulf + \h}{d, - + h,h,_ih. (3.3b) 

and 

n 

h ■■= + + h,(^% - + -I- ^ h~d,, (3.4a) 

6^=u, •< + (/i,_i + i/i,)C, (3.4b) 

which are the mean and a components of the ith layer Bernoulli head. The above re- 
sult follows upon realizing that ^Yl!]= \ hjdjdthj^i — ^thj X]fc=j+i ^k^k = 0, and is 
largely facilitated by rewriting (|2.7t l.c') in the form 

dtUi + fi^uf + hiz X (5iUj + uf ) + Vbi = R^, (3.5a) 

9t< + Ki X (gf u, + q,n1) + Vh^ = (3.5b) 

Here, 

:= i/iriy • h,xi1 (3.6) 
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is the vertical average of the ith cr-vertical velocity; 



:= hr^ (/ + V • u, X z) , hr^V • < x z (3.7) 

are the mean and a components of the ith-layer tr-potential vorticity, which are not 
Lagrangian constants of (|2.7|l : and 

R, Vi?, + ^h,Vid, - K), (3.8a) 

:= {h-i + \h,)Vr^ - ^h,V^„ (3.8b) 

which are rotational forces that arise as a consequence of the buoyancy inhomogeneities 
within each layer {V^, 7^0^ ^^1)- 

In turn, the local conservation law for the sum of the zonal momenta within each layer 
is given by 

dt + V • ^ Ff + d.,ho hj^J = 0, (3.9a) 

J j j 

where 

i-l 

Ff := M,u, + \h,u1xi1 + \-ihl{{>, - + -ih,+iih+i Y (3.9b) 

i=i 

with Ui denoting the zonal component of and x the unit vector in the same directionf . 
The above result follows upon multiplying by jhi the zonal component of (|2.7H ). 

^tu^ + u, • Vu, + ^hr^V ■ hulnl - (/ + tu,) v, + 7"^a;^» = 0, (3.10) 

and realizing that Y^j^i^ j^^Ahj-i - ha) + Y.j^jdxYJl^j+ihk'&k = d^ij^j f^j+i'^j+i 
hk). At this point it is crucial to specify whether the geometry is flat or spherical. 
On the sphere, Va = (j~^dxa, dyOj , for any scalar a(x), and V • a = j^^[dxa + dy (7^)], 
for any vector a — (a, b), where a; = (A — Aq) cos 9 R and y = [6 — 9q)R are, respectively, 
rescaled geographic longitude and latitude on the surface of the Earth whose mean radius 
is R\ and 7 — cos 6*0 cos and r — R^^ta.nd = —^~^d'y /dy are coefficients that char- 
acterize the geometry of the space (the arclength element square and area element are 
dx^ = 7^da;^ + dy^ and d^x = jdxdy, respectively). The ith zonal momentum (angular 
momentum around the Earth's axis) is then given by 

M, = h.l-yu, - rji?(cosi?o - 7 cos i?)], (3.11) 

where il is the Earth's angular rotation rate. In the classical (3 plane, 7=1 and t = 
so that all the operators are Cartesian and Mi = hi{ui — foy — ^Py'^). However, the 
geometry in a consistent /3 plane cannot b e Cartesian : instead 7 = 1 — Toy, r = T0/7, 
and Mi — hil-fUi — foy — i/3(l — R'^r^y'^] l)Ripalll99'7|) . Conservation of the total zonal 
momentum (sum over all layers) in a closed domain requires in addition, in all cases, 
that both the topography and coasts to be zonally symmetric. 

f A term —^dxih'&a) is missing on the right hand side of (4.6) in R95. 
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3.2. Kelvin circulation theorem(s) 

Letting Ct{ui) be a material curve that moves in the ith layer with the vertically averaged 
velocity u,;. Then taking into account H2.7b ). from (??) it follows that 



f dx- (u, -l-u/) = i dx- (R, + X z + /i,0, (3.12a) 

(f> dx-uf=i dx • (Rf + /i,9,uf X z) , (3.12b) 



d 

di 

d_ 

dt 



where V • u/ x z := /. These are the Kelvin circulation thcorem(s) for the multilayer 
model of this paper. The situation is different than in the fully three-dimensional model 
for which fdx • (v|^ + Uj), where v|^ is the (three-dimensional) velocity on a given 
isopycnal surface, is conserved. In the present case, Ct{ui) cannot remain on an isopycnal 
surface, which moves with u|^ at all times. If instead of Ct{ui) a rigid boundary were 
considered on which n x Vi?i = = n x Vi?°^, where fi is the outer normal to the 
boundary, then the circulation of u°' along that boundary would be preserved. Even 
though preservation of such a boundary condition is guaranteed by the dynamics, it is 
too restrictive to be expected. In the slab model, however, fi x Vi?i = on a rigid 
boundary, which is preserved by the dynamics, is necessary for the circulation along that 
boundary to be an (explicit) integral of motion, and for the Hamiltonian structure of the 
model not to be spoiled. 



3.3. Hamiltonian structure 

The (inviscid unforced) fully three-dimensional primitive equations model, just as the 
Euler equations of fluid mechanics, possesses what is called in geometric mechanics a 
"generalized or singular Hamiltonian structure." A good sign of the validity of any ap- 
proximate model derived from it is the preservation of such a structure. This topic is 
treated here; readers not interested in geometric mechanics may safely skip this section. 

Upon grouping (p = (??i, , /ii, Ui, uj, • • • ,■!?„, i?;^ , /i„ , u„ , <) and defining H[ip] := 
J2j Jd d^xE'j, where the model equations 12. 7|) can be cast in the form 

dt^ = {v,n}. (3.13) 

Here, {A, B} :— d^x {SA/Sip)J{SB/Sip) for any admissible functionals A[ip] and B[(p], 
where J is a skew-adjoint 7x7 block-diagonal matrix operator J — 0^ J^^.^ + with 
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gfzx 







qfzx 



(3.14a) 
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-3/i-iVi9, 
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V(/i,-i< 








(3.14b) 



Operators 1)3.14(1 are the same ones given by R95 for the one-layer model but particularized 
for ith layer. Thus an admissible functional A[(f] must be such that {6A/6ui) • fi = = 
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(SA/Suf) ■ n on dD. The functional derivatives of Ti are given by 
6n , ,~ 1, ^ -57^ SH , r , _ 6H 

OU,- ^ 017, 6 Ohi OUj ouf 



That system (12 .Tf) can be put in the form (|3.13l) suggests that the layered model of this 
paper posses a singular Hamiltonian structure on the phase space given by the Cartesian 
product of the n manifolds whose coordinates are the 7 layer variables. The functional H 
and the operator {•, •} being the Hamiltonian and Poisson bracket, respectively. (Proof of 
the Jacobi identity satisfied by the Poisson bracket — not given in R95 — will be published 
elsewhere.) This conveys important properties like the direct linking of the conservation 
laws for energ y, ?i, and zonal momentum, A4, with symmetries via Noether's theorem 
(e.g. Shcphcrdl ll99q) : notice that Ti and —M are the generators of t- and x-translations 
because of (|3.13|) and d^if = y?}, respectively. In this framework the ith layer inte- 
grals of volume, , mass, C}, and buoyancy variance, Cf, are invariant Casimirs, which 
satisfy {Cf,J^} = yT[ip] and thus are not related to (exphcit) symmetries; notice that 
Cf does not generate any transformation since {(p,Cf} = 0. 

The existence of a Lie-Poisson Hamiltonian structure, which would imply an analogous 
Euler-Poincare structure, has not been established for the present model. The Euler- 
Poincare structure is based on a Lagrangian functional a nd Hamilton's pri nciple, rather 
than a Hamiltonian functional and a Poisson bracket I Holm et a/.ll200^ . The use of 



constrained dynamics 


re.g.lSa,]moJl983HHolmlll99fil: 


Beron- Veral2003^ and approximate 


dynamics l|Holm et al\ 


20021). An interesting research avenue is the exploration of whether 



approximations directly made in the fuUy-three dimensional Lagrangian can lead to the 
approximate model of the present paper. The study of these topics are reserved for the 
future. 

3.4. Arnold stability 

In R95 it was shown that for one layer a state of rest (or a steady state with at most a 
uniform zonal current) can be shown to be formally stable using Arnold's (1965, 1966) 
method if and only if (|2.6I) is satisfied, i.e. if and only if the buoyancy is everywhere 
positive and increases (resp., decreases) with depth within a layer with the rigid bottom 
(resp., rigid lid). Arnold's method for proving the stability of steady solution of a system 
consists in searching for conditions that guarantee the sign-definiteness of a general invari- 
ant which is quadratic to the lowest- order in the deviat i on from that state; the resul ting 
conditions being only sufficient (e.g. iHolm et a^jri985l iMcIntvre fc Shepherdlll987j) . In 
the present model, however, Arnold's method fails to provide stability conditions even 
for a state of rest and with no topography {ho = 0). The lowest-order contribution to 
that quadratic invariant, which can be called a "free energy" because it corresponds to 
a state of rest, 



+ {gjShj + Hjddj)dhj_i, (3.16) 
cannot be proved sign-definite. Here, 

I d'x(-), (3.17) 
3 , -Id, 
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and Hi, gi and iVi, which are the ith-layer unperturbed depth, verticaUy averaged buoy- 
ancy, and Brunt-Vaisala frequency, respectively. Similarly, a state of rest in the multilayer 
slab model cannot be proved formally stable using Arnold's method. Surprisingly, it is 
possible to prove the stability of a steady state with a uniform zonal current in that 
model. But the condition of stability is not one of "static" stability like H2.6|l as in the 
one-layer version of Ripa's model. Contrarily, it is one of "baroclinic" stability since 
a uniform current in the multilayer slab model has an implicit vertic al shea r th rough 
the th ermal-wind balance. These results can all be inferred from iRipal l|l993|) and lRipal 
lll99fil) . 

It is worthwhile to mention, however, that there is at least a system, which has one 
Ripa-like layer and n — 1 homogeneous layers, for which a state of rest can be proved 
formally stable. For instance, assuming the uppermost layer to be Ripa-like, the corre- 
sponding free energy takes the form 



+ fe - 9i+i){5hj f - \NlH^{5h^)\ (3.18) 

where a := n (resp., a := 1) for the rigid-bottom (resp., rigid-lid) configuration, and 
Hi, gi, and Ni are all constants. The above free energy is positive-definite if and only if 
1)2. 6|l if fulfilled. [The homogeneous layers model has an infinite set of invariants which 
are given by Jj hjF(qj) Vi^(-); these include the volume integral, which is the only one 
needed to obtain the above result.] When all layers are homogeneous the same result is 
obtained. When one slab layer is included, however, the free energy cannot be shown of 
one sign. 

The fact that it is not possible to prove formal stability for a steady state in the 
multilayer model of this paper does not mean that state is unstable; it simply means 
that Arnold's method is not useful to prove it. One reason for the failure of Arnold's 
method is the lack of vorticity-related invariants. 

3.5. Waves and instabilities 

In R95 it was shown that a model with one layer supports two vertical normal modes. 
The present n-layer model can support 3n — 1 in account for the extra n — I baroclinic 
modes contributed by the n layers. A model with n homogeneous or slab layers only 
supports 2n vertical normal modes. Figure shows, in particular, the phase speed of 
long gravity waves (c) for the barotropic and first baroclinic mode as a function of the 
stratification in the reference state, as predicted by the (exact) fully three-dimensional 
model, 

NrH N,H, S 

2c c 1 — 5 

(cf. e.g. lGilll|l982|) . and as predicted by the present model with one and two layers, and 
also by a model with two homogeneous layers. Here, 

S \N^H,/g, (3.20) 

is a nondimensional measure of the stratification within the active layer, where H^, g,-, and 
TVr denote the constant layer thickness, vertically averaged buoyancy, and Brunt-Vaisala 
frequency, respectively, in the motionless reference state (R95, Beron-Vera & Ripa, 1997). 
The reference state has a flat bottom/lid {Hq = 0) and a single, continuously stratified 
active layer with linear vertical stratification. The buoyancy varies from ^^(l =p 5*), at the 
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Figure 3. Phase speed of a long gravity wave for the barotropic (cq) and first barochnic (ci) 
modes as a function of the stratification of the reference state, as predicted by the fuUy-three 
dimensional model (solid lines), a one- (dot-dashed lines) and two-layer (dashed) versions of 
the present multilayer model, and a two-homogenous-layer model (dotted lines). In this figure, 
Hr, Qr, and Nr are the constant reference depth, mean buoyancy, and Brunt-Vaisala frequency, 
respectively. 



top of the active layer, to 5r(l ± 5"), at the base of the active layer. As a consequence of 
condition (|2.6II . Q < S < 1. Notice that this condition does not hmit the size of the density 
gradient within the active layer: the buoyancy change within the layer, N^H^ = 2Sgr, 
can be made as large as desired while keeping fixed the buoyancy at the interface with 
the inert layer g^-^l — S). The present layer model reference state parameters are chosen 
as 

2i-V 



Hi = — , ft 



1+1 



N. = ^, (3.21) 



with n = 1,2. The ith layer reference depth Hi and layer buoyancy gi for the model with 
two homogeneous layers are chosen as in (|3.21|l . [Other choices for the two-homogeneous- 
layer model parameters, which can be made to match the three-dimensional model 's 
phase speed at a fixed valu e of S, are of course possible (e.g. lOlascoaea fc RipalHoogfl .] 
For this model it is (cf. e.g. lGil]|ll982(l 



c2 ^ g^H, S S /2 

The long-gravity-wave phase speeds of the model of this paper are computed numerically 
by setting / = in (|2.7(l and considering infinitesimal normal-mode perturbations to 
l|3.21|l . The left-hand-side panel of figure 01 shows that a one-layer version of the model 
generalized here (dot-dashed line) is enough to represent the exact (solid line) barotropic 
mode phase speed (cq) for all values of the stratification parameter S (R95). The result 
from a two-layer version (dashed line) cannot be distinguished from the one-layer version 
and exact results in this plot. The two-homogeneous-layer model (dotted line) gives the 
exact phase speed at very weak stratification {S — > 0) but overestimates it somewhat as 
the stratification increases. The most striking differences appear in the representation of 
the baroclinic mode phase speed. The right-hand-side panel of figure 13 shows that the 
one-layer version of the model of this paper provides only a fair representation of the 
phase speed of the gravest baroclinic mode (ci). The two-layer version, instead, provides 
an excellent approximation for any size of 5*. The two-homogeneous-layer model phase 
speed, in turn, overestimates the exact model phase speed for all S (the absolute value 
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of difference being of the same order as that produced by the one-layer version of the 
present model). The reason for the improved performance of the two-layer version of 
present model compared to the one-layer version and the two-homogeneous-layer model 
is that with one linear function in each layer is possible to better represent the exact 
vertical structure eigenfunctions than with only one linear function and one number per 
layer, respectively. A considerable amount of "physics" is gained by adding to the ve- 
locity and buoyancy fields a vertical sh ear within each layer following the simple recipe 
proposed by R95. In iBeron-Vera et "dl (1200.1) we extend and complete the above calcu- 
lation by considering more complicated reference vertical stratifications and computing 
the vertical structure eigenfunctions. We also show that the present model supports the 
usual midlatitude and equatorial gravity and vortical waves (Poincare, Kelvin, Rossby, 
Yanai, etc.) in each of the vertical normal modes. 

Let me close this section by making some remar ks on baroclinic in stability results. 
To study baroclinic instability in the upper ocean, iBeron-Vera fc Rioa (,1997.1 consid- 
ered a reduced-gravity approximation. In that work it was shown that a reduced-gravity 
one-layer version of Ripa's model gives, in the limit of weak internal stratification and 
sufficiently strong velocity's vertical shear, the exact solution for normal-mode pertur- 
bations whose scales are of the order of the external (equivalent barotropic) deformation 
radius of the system. For perturbations of the order of th e internal (first baroclinic) 
deformation radius, i.e. the classical Eady "waves" (cf. e.g. IPedloskvlll987l) . the model 
resulted less successful as the exact eigensolutions have an exponential trapping behav- 
ior which can be only poorly represented if the velocity and buoyancy vertical profiles 
are restricted to a linear function of depth. The very important dynamical aspect of 
the presence of a high-wavenumber cutoff of instability was nonetheless captured by the 
one-layer model, though overestimated in about 40 %. By virtue of the above results, I 
expect this overestinia tion to be substantially r educed if two active layers are considered. 
This is investigated in iBeron- Vera et all l)2003(l where we further test the present model 
performance in upper-ocean baroclinic instability including ageostrophic eff ects. The lat- 
ter is done by relaxing the rigid bottom boundary constraint in the classical lstond l)l966() 



In R95 forcing (wind stress, interfacial drag, and buoyancy/heat input) was introduced 
in the one-layer version of the present multilayer model equations in a way that was 
compatible with the conservation laws of energy, momentum, and mass/heat content. 
The same approach is adopted here to include, in addition, freshwater fluxes through the 
surface in accordance with the conservation law of salt content. The possibility for the 
exchange of fluid across the other interfaces is also considered. 

Let r(x, t) be a wind stress acting at the surface of the ocean {Pn+i = must be 
the setting in the rigid-bottom configuration and typically /ig = in the rigid- lid one). 
Assume further that there is a friction force acting at the interface between contiguous 
layers. Introduction of these forces in Newton's equations (|2.7t l.e) in the form 



implies that the work done by the wind stress is proportional to the velocity at the top 
of the uppermost layer, Uc =F Uq, and that one done by the friction force in the ith layer 



model. 



3.6. Forcing 



dtUi H = SiaT/ha - ri{Ui ± u^), 

9tuf H = T55.,aT/ha + 3r^(ui ± u^. 



(3.23a) 
(3.23b) 
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is proportional to the velocity at the base of that layer, u; ± uf . Namely 

j 3 

dtJ2M,+--- = T-i,-J2 rjh.iu, ± uj) • X. (3.24b) 

j j 

In the above equations, 5ij is the Kroenecker delta and is a friction coefficient that 
can be taken as a constant or as some function of hi and |ui ± uf |. [Recall that a := n 
(resp., a :— 1) for the rigid-bottom (resp., rigid-lid) configuration.] 

Let now r(x, t) be a buoyancy input through the surface and introduce it in the 
buoyancy equations (|2.7b .bl in the form 

dtd, + ---^ S,^r/h^, (3.25a) 
atz?^ + • • • - r]S,c.r/h^, (3.25b) 

where rj is any constant. Consider, in addition, the possibility of fluid crossing the interface 
between consecutive layers; then the volume conservation equation (|2.7b ') can be rewritten 
as 

dth, + -- - ^w^ -wl (3.25c) 

Here, the quantities (x, t) and (x, t) are volume fluxes per unit area through the top 
and base of the ith layer, respectively. The set (|3.25|) . for any value of r/, is compatible 
with the mass conservation equation 

^tih^^^) + ■■■ = S,^r + - wl). (3.26) 

At the surface wl^{x, t) = P(x, t)—E{x, t), which represents the imbalance of precipitation 
minus evaporation; away from the surface some param etrization mu st be adopted. In 
models with slab layers it is commonly set (e.g. Mc Crcarv et a^jll99l|) 

-wl^l ^-^y^'^^^^f^^ ^ ^-1' (3.27) 

[ if > i7<Li. 

Here, H° and are constants that with units of length and time, respectively, that 
characterize the "entrainment" process. In the present case, an algorithm may be design 
such that condition (|2.6(l is fulfilled at all times. This would allow for the inclusion of 
mixing processes in a novel and rather physical way. On the other hand, mixing events 
are known to occur in localized regions of the ocean, which can be represented by regions 
where -I- ^9°^^ < di — instantaneously. This subject deserves to be studied in 
detail. 

Let finally assume a linear state equation, i.e. 1?^ = gariTi — Tn+i) — gas{Si — Sn+i). 
Here, ax and as are the thermal expansion and salt contraction coefficients, respec- 
tively; T,{x, cr, t) = Ti(x, t) + (jT[{yL, t) and S',(x, cr, t) = ^^(x, t) + crS^{yi, t) are the ith 
layer temperature and salinity, respectively; and r„+i and Sn+i are the inactive layer 
(constant) temperature and salinity, respectively. Let also write the buoyancy input as 

F - gaT{.P,Cp)-^Q + gagS^iP ~ E). (3.28) 

Here, Cp is the specific heat at constant pressure and Q{'x.,t) is the heat input through 
the surface. Equation H3.26|l can then be split into a heat and salt content conservation 
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equations, namely 

dtih,f,) + • • • - S,^ip^Cp)-^Q + T,;K^ - wl), (3.29a) 
dt{h,S^) + • • • - 6^ME -P) + S,iw'^ - wl). (3.29b) 

If fluid across the surfa ce is allowed only, the choice 13.28|l enforces, on one hand (e.g. 
iBeron-Vera et al\\l99^ 

^ / h,S, ^ 0, (3.30a) 



dt 

- J 

and, on the other iBeron- Vera fc Ripal ll200nt) 

^ (T) = d^x (p,Cp)-ig + (T„ - (T))(P - E), (3.30b) 

where ^ :— ^. hj = d'^x h is the total volume and (T) := J^. hjTj is the aver- 
age temperature in V. Notice that H3.30b ). unlike the equation satisfied b y h jTj, is 
independent of the choice of the origin of the temperature scale fcf. IWarrerJll999|) . 



4. Concluding remarks 

This paper describes a multilayer extension of Ripa's (1995) primitive equations one- 
layer ocean model with thermodynamics. Inside each layer the velocity and buoyancy 
fields can vary not only arbitrarily in the horizontal position and time, but also linearly 
with depth. In the absence of external forcing and dissipation, the model conserves vol- 
ume, mass, buoyancy variance, energy, and zonal momentum for a zonally symmetric do- 
mains and topographies. Furthermore, as the fully three-dimensional primitive equations 
from which it is derived, the present model possesses a singular Hamiltonian structure. 
Unlike the multilayer slab models, i.e. those for which the velocity and buoyancy fields 
are depth-independent within each layer, the model generalized here is able to represent 
the thermal-wind balance explicitly at low frequency inside each layer. In this sense, the 
model of this paper has better "physics" than a model with slab-like layers. For the 
same number of layers, n, say, on the other hand, the model of this paper can sustain 
2n — 1 more vertical normal modes as homogeneous layers models, which, in addition, 
are not able to incorporate thermodynamic processes. A model with only two layers has 
been shown here to provide an excellent representation of the exact first baroclinic mode 
phase speed. In this other sense, the present model has more "physics" than a model 
with homogeneous layers. The present generalization enriches Ripa's single-layer model 
by providing it enough flexibility to approach problems for which a single-layer structure 
is too idealized. Configurations with a small number of layers are particularly useful for 
the insight they provide into physical processes. Configurations with more layers can be 
used as the basis for an accurate numerical circulation model. 



This paper has benefited from conversations with Josefina Olascoaga, Javier Zavala- 
Garay, Mohamed Iskandarani, Mike Brown, and Rigo Garcia. 
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